We present a rice-wheat systen spatial decisionb support tool which relies on data from the matched landscape crop assessment surveys in eastern India and a multivariate geoadditive Bayesian model to predict entry points for system optimization in the area of interest.

Code 
library(sp)
library(mgcv)
library(bamlss)

# Multivariate geoadditive model
# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")
# library(mvnchol)
library(BayesX)
library(R2BayesX)
library(sf)
library(spdep)
library(rio)
Irrig_Rev_rice_wheat <- import("data/Irrig_Rev_rice_wheat.csv")

library(janitor)
# Irrig_Rev_rice_wheat <- clean_names(Irrig_Rev_rice_wheat)
library(tidyr)
#Irrig_Rev_rice_wheat <- Irrig_Rev_rice_wheat %>% drop_na()

# library(lubridate)
# library(anytime)
# Irrig_Rev_rice_wheat$januaryfirst2017 <- ymd("2017-01-01")

# Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day <- Irrig_Rev_rice_wheat$sowdate_fmt - Irrig_Rev_rice_wheat$januaryfirst2017

# Irrig_Rev_rice_wheat$sowdate_fmt_rice_day <- Irrig_Rev_rice_wheat$sowdate_fmt_rice - Irrig_Rev_rice_wheat$januaryfirst2017

Irrig_Rev_rice_wheat$harvest_day_rice <- Irrig_Rev_rice_wheat$l_crop_duration_days_rice + Irrig_Rev_rice_wheat$sowdate_fmt_rice_day

# Irrig_Rev_rice_wheat$harvest_day_wheat <- Irrig_Rev_rice_wheat$l_crop_duration_days + Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day


shpname <- file.path(getwd(), "shp", "India_aoi_sf_sp")
India_aoi_sp_bnd <- BayesX::shp2bnd(shpname = shpname, regionnames = "District", check.is.in = F)
Reading map ... finished
Note: map consists originally of 50 polygons
Note: map consists of 47 regions
Code 
f_rice_yield_MRF <- list(
    b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sigma ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)



K <- neighbormatrix(India_aoi_sp_bnd)
head(K)
           Araria Arwal Aurangabad Banka Begusarai Bhagalpur Arah Buxar
Araria          4     0          0     0         0         0    0     0
Arwal           0     6         -1     0         0         0   -1     0
Aurangabad      0    -1          3     0         0         0    0     0
Banka           0     0          0     3         0        -1    0     0
Begusarai       0     0          0     0         5         0    0     0
Bhagalpur       0     0          0    -1         0         6    0     0
           Darbhanga Gaya Gopalganj Jamui Jehanabad Kaimur Katihar Khagaria
Araria             0    0         0     0         0      0       0        0
Arwal              0   -1         0     0        -1      0       0        0
Aurangabad         0   -1         0     0         0      0       0        0
Banka              0    0         0    -1         0      0       0        0
Begusarai          0    0         0     0         0      0       0       -1
Bhagalpur          0    0         0     0         0      0      -1       -1
           Kishanganj Lakhisarai Madhepura Madhubani Munger Muzaffarpur Nalanda
Araria             -1          0        -1         0      0           0       0
Arwal               0          0         0         0      0           0       0
Aurangabad          0          0         0         0      0           0       0
Banka               0          0         0         0     -1           0       0
Begusarai           0         -1         0         0     -1           0       0
Bhagalpur           0          0        -1         0     -1           0       0
           Nawada WestChamparan Patna EastChamparan Purnia Rohtas Saharsa
Araria          0             0     0             0     -1      0       0
Arwal           0             0    -1             0      0     -1       0
Aurangabad      0             0     0             0      0     -1       0
Banka           0             0     0             0      0      0       0
Begusarai       0             0    -1             0      0      0       0
Bhagalpur       0             0     0             0     -1      0       0
           Samastipur Saran Sheikhpura Sheohar Sitamarhi Siwan Supaul Vaishali
Araria              0     0          0       0         0     0     -1        0
Arwal               0     0          0       0         0     0      0        0
Aurangabad          0     0          0       0         0     0      0        0
Banka               0     0          0       0         0     0      0        0
Begusarai          -1     0          0       0         0     0      0        0
Bhagalpur           0     0          0       0         0     0      0        0
           Balia Chandauli Deoria Gazipur Gorakhpur Kushinagar Maharajganj Mau
Araria         0         0      0       0         0          0           0   0
Arwal          0         0      0       0         0          0           0   0
Aurangabad     0         0      0       0         0          0           0   0
Banka          0         0      0       0         0          0           0   0
Begusarai      0         0      0       0         0          0           0   0
Bhagalpur      0         0      0       0         0          0           0   0
           Siddharthnagar
Araria                  0
Arwal                   0
Aurangabad              0
Banka                   0
Begusarai               0
Bhagalpur               0
Code 
## Also need to transform to factor for
## setting up the MRF smooth.
Irrig_Rev_rice_wheat$District <- as.factor(Irrig_Rev_rice_wheat$a_q103_district)

## Now note that not all regions are observed,
## therefore we need to remove those regions
## from the penalty matrix
rn <- rownames(K)
lv <- levels(Irrig_Rev_rice_wheat$District)
i <- rn %in% lv
K <- K[i, i]

set.seed(321)
library(bamlss)
b_rice_yield_MRF <- bamlss(f_rice_yield_MRF, data = Irrig_Rev_rice_wheat, family = "gaussian")
AICc 12433.92 logPost -6216.75 logLik -6063.08 edf 148.79 eps 2.3399 iteration   1
AICc 12239.92 logPost -5910.29 logLik -5967.69 edf 147.28 eps 0.4613 iteration   2
AICc 12211.69 logPost -5884.01 logLik -5956.93 edf 144.14 eps 0.1461 iteration   3
AICc 12205.11 logPost -5882.70 logLik -5953.35 edf 144.42 eps 2.9473 iteration   4
AICc 12203.12 logPost -5882.46 logLik -5952.28 edf 144.48 eps 0.0429 iteration   5
AICc 12202.50 logPost -5882.43 logLik -5951.95 edf 144.50 eps 0.0113 iteration   6
AICc 12202.32 logPost -5882.45 logLik -5951.84 edf 144.52 eps 0.0088 iteration   7
AICc 12202.28 logPost -5882.44 logLik -5951.82 edf 144.52 eps 0.0010 iteration   8
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0002 iteration   9
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0001 iteration  10
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0001 iteration  11
AICc 12202.26 logPost -5882.44 logLik -5951.80 edf 144.53 eps 0.0000 iteration  12
AICc 12202.26 logPost -5882.44 logLik -5951.80 edf 144.53 eps 0.0000 iteration  12
elapsed time:  7.82sec
Starting the sampler...

|                    |   0%  1.19min
|*                   |   5%  1.21min  3.83sec
|**                  |  10%  1.15min  7.66sec
|***                 |  15%  1.10min 11.66sec
|****                |  20%  1.05min 15.77sec
|*****               |  25%  1.01min 20.13sec
|******              |  30% 57.10sec 24.47sec
|*******             |  35% 53.63sec 28.88sec
|********            |  40% 49.71sec 33.14sec
|*********           |  45% 45.76sec 37.44sec
|**********          |  50% 41.75sec 41.75sec
|***********         |  55% 37.70sec 46.08sec
|************        |  60% 33.66sec 50.49sec
|*************       |  65% 29.48sec 54.75sec
|**************      |  70% 25.32sec 59.07sec
|***************     |  75% 21.13sec  1.06min
|****************    |  80% 16.94sec  1.13min
|*****************   |  85% 12.73sec  1.20min
|******************  |  90%  8.48sec  1.27min
|******************* |  95%  4.23sec  1.34min
|********************| 100%  0.00sec  1.41min
Code 
## First, note that we have the structured id = 'mrf1' and unstructured
## spatial effect id = 're2', also indicated in the model summary
summary(b_rice_yield_MRF)

Call:
bamlss(formula = f_rice_yield_MRF, family = "gaussian", data = Irrig_Rev_rice_wheat)
---
Family: gaussian 
Link function: mu = identity, sigma = log
*---
Formula mu:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + 
    s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + 
    s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              1.90903 0.53256 1.64495 3.75492      0.013
rice_duration_class_long 0.07427 0.01184 0.07438 0.13788      0.075
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(sowdate_fmt_rice_day).tau21     5.510e-01 1.460e-02 2.985e-01 2.446e+00
s(sowdate_fmt_rice_day).alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(sowdate_fmt_rice_day).edf       3.784e+00 1.625e+00 3.700e+00 5.909e+00
s(g_q5305_irrig_times_rice).tau21 2.324e+00 2.653e-01 1.653e+00 8.760e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf   6.482e+00 4.641e+00 6.530e+00 7.994e+00
s(nperha_rice).tau21              4.816e-01 1.694e-04 7.272e-02 3.602e+00
s(nperha_rice).alpha              1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf                2.866e+00 1.010e+00 2.509e+00 6.675e+00
s(p2o5perha_rice).tau21           9.498e-02 9.163e-05 7.279e-03 8.451e-01
s(p2o5perha_rice).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf             1.889e+00 1.007e+00 1.411e+00 4.674e+00
s(District,id='mrf1').tau21       4.290e-02 7.425e-05 1.089e-02 2.036e-01
s(District,id='mrf1').alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf         2.274e+01 1.975e+00 2.845e+01 3.410e+01
s(District,id='re2').tau21        5.664e+00 2.034e-01 5.310e+00 1.506e+01
s(District,id='re2').alpha        1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf          3.560e+01 3.391e+01 3.581e+01 3.594e+01
                                  parameters
s(sowdate_fmt_rice_day).tau21          0.182
s(sowdate_fmt_rice_day).alpha             NA
s(sowdate_fmt_rice_day).edf            3.291
s(g_q5305_irrig_times_rice).tau21      2.895
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        7.106
s(nperha_rice).tau21                   1.398
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     5.489
s(p2o5perha_rice).tau21                0.018
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  1.787
s(District,id='mrf1').tau21            0.029
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             31.761
s(District,id='re2').tau21            12.212
s(District,id='re2').alpha                NA
s(District,id='re2').edf              35.917
---
Formula sigma:
---
sigma ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + 
    s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + 
    s(District, bs = "mrf", xt = list(penalty = K)) + s(District, 
    bs = "re")
-
Parametric coefficients:
                              Mean      2.5%       50%     97.5% parameters
(Intercept)              -0.119970 -0.179496 -0.119649 -0.063475     -0.114
rice_duration_class_long  0.061685  0.009405  0.061474  0.115939      0.060
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9860 0.8949 1.0000     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(sowdate_fmt_rice_day).tau21     2.255e-02 9.588e-05 4.358e-03 1.617e-01
s(sowdate_fmt_rice_day).alpha     9.779e-01 8.326e-01 1.000e+00 1.000e+00
s(sowdate_fmt_rice_day).edf       1.623e+00 1.009e+00 1.340e+00 3.522e+00
s(g_q5305_irrig_times_rice).tau21 2.200e-01 2.578e-04 1.086e-01 1.030e+00
s(g_q5305_irrig_times_rice).alpha 9.450e-01 6.660e-01 9.945e-01 1.000e+00
s(g_q5305_irrig_times_rice).edf   4.016e+00 1.123e+00 4.209e+00 6.532e+00
s(nperha_rice).tau21              1.741e-01 4.189e-04 3.295e-02 1.329e+00
s(nperha_rice).alpha              9.529e-01 6.319e-01 9.979e-01 1.000e+00
s(nperha_rice).edf                2.684e+00 1.043e+00 2.352e+00 6.037e+00
s(p2o5perha_rice).tau21           3.062e-01 1.414e-04 1.220e-02 2.836e+00
s(p2o5perha_rice).alpha           9.579e-01 6.586e-01 9.982e-01 1.000e+00
s(p2o5perha_rice).edf             2.461e+00 1.017e+00 1.808e+00 6.536e+00
s(District,id='mrf1').tau21       2.374e-03 7.354e-05 1.646e-03 7.531e-03
s(District,id='mrf1').alpha       8.160e-01 2.828e-01 9.140e-01 1.000e+00
s(District,id='mrf1').edf         1.912e+01 2.882e+00 2.035e+01 2.921e+01
s(District,id='re2').tau21        2.233e-02 5.515e-03 2.132e-02 4.458e-02
s(District,id='re2').alpha        7.382e-01 1.551e-01 8.258e-01 1.000e+00
s(District,id='re2').edf          2.842e+01 1.981e+01 2.919e+01 3.199e+01
                                  parameters
s(sowdate_fmt_rice_day).tau21          0.014
s(sowdate_fmt_rice_day).alpha             NA
s(sowdate_fmt_rice_day).edf            1.802
s(g_q5305_irrig_times_rice).tau21      0.191
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        4.768
s(nperha_rice).tau21                   0.091
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     3.139
s(p2o5perha_rice).tau21                2.253
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  6.288
s(District,id='mrf1').tau21            0.007
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             28.648
s(District,id='re2').tau21             0.002
s(District,id='re2').alpha                NA
s(District,id='re2').edf              10.540
---
Sampler summary:
-
DIC = 12123.05 logLik = -6011.865 pd = 99.3237
runtime = 85.53
---
Optimizer summary:
-
AICc = 12202.27 edf = 144.5353 logLik = -5951.809
logPost = -5882.444 nobs = 4537 runtime = 7.82
Code 
# Plot the nonlinear effect
plot(b_rice_yield_MRF, model = "mu", term = "s(sowdate_fmt_rice_day)")
Code 
plot(b_rice_yield_MRF, model = "mu", term = "s(g_q5305_irrig_times_rice)")
Code 
plot(b_rice_yield_MRF, model = "mu", term = "s(nperha_rice)")
Code 
plot(b_rice_yield_MRF, model = "mu", term = "s(p2o5perha_rice)")
Code 
f_wheat_yield_MRF <- list(
    l_ton_per_hectare ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sigma ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)
set.seed(321)
b_wheat_yield_MRF <- bamlss(f_wheat_yield_MRF, data = Irrig_Rev_rice_wheat, family = "gaussian")
AICc 8958.541 logPost -4419.70 logLik -4335.18 edf 139.61 eps 0.3246 iteration   1
AICc 8445.144 logPost -4096.84 logLik -4071.93 edf 145.76 eps 0.1431 iteration   2
AICc 8344.698 logPost -4042.08 logLik -4030.62 edf 137.40 eps 0.0345 iteration   3
AICc 8319.455 logPost -4032.32 logLik -4022.32 edf 133.33 eps 0.0132 iteration   4
AICc 8313.001 logPost -4030.95 logLik -4020.60 edf 131.91 eps 0.0050 iteration   5
AICc 8312.046 logPost -4029.90 logLik -4020.23 edf 131.81 eps 0.0019 iteration   6
AICc 8311.784 logPost -4029.54 logLik -4020.15 edf 131.76 eps 0.0010 iteration   7
AICc 8311.335 logPost -4029.64 logLik -4020.12 edf 131.57 eps 0.0005 iteration   8
AICc 8311.297 logPost -4029.53 logLik -4020.11 edf 131.57 eps 0.0003 iteration   9
AICc 8311.268 logPost -4029.43 logLik -4020.10 edf 131.56 eps 0.0002 iteration  10
AICc 8311.245 logPost -4029.35 logLik -4020.09 edf 131.56 eps 0.0002 iteration  11
AICc 8311.226 logPost -4029.28 logLik -4020.09 edf 131.56 eps 0.0002 iteration  12
AICc 8311.211 logPost -4029.22 logLik -4020.08 edf 131.56 eps 0.0001 iteration  13
AICc 8311.197 logPost -4029.16 logLik -4020.08 edf 131.55 eps 0.0001 iteration  14
AICc 8311.185 logPost -4029.11 logLik -4020.07 edf 131.55 eps 0.0001 iteration  15
AICc 8311.175 logPost -4029.06 logLik -4020.07 edf 131.55 eps 0.0001 iteration  16
AICc 8311.175 logPost -4029.06 logLik -4020.07 edf 131.55 eps 0.0001 iteration  16
elapsed time:  5.11sec
Starting the sampler...

|                    |   0% 52.36sec
|*                   |   5% 53.77sec  2.83sec
|**                  |  10% 50.58sec  5.62sec
|***                 |  15% 48.51sec  8.56sec
|****                |  20% 46.48sec 11.62sec
|*****               |  25% 44.49sec 14.83sec
|******              |  30% 42.54sec 18.23sec
|*******             |  35% 40.10sec 21.59sec
|********            |  40% 38.51sec 25.67sec
|*********           |  45% 37.64sec 30.80sec
|**********          |  50% 36.25sec 36.25sec
|***********         |  55% 33.61sec 41.08sec
|************        |  60% 30.01sec 45.02sec
|*************       |  65% 26.28sec 48.80sec
|**************      |  70% 22.53sec 52.56sec
|***************     |  75% 18.74sec 56.23sec
|****************    |  80% 15.02sec  1.00min
|*****************   |  85% 11.28sec  1.07min
|******************  |  90%  7.54sec  1.13min
|******************* |  95%  3.78sec  1.20min
|********************| 100%  0.00sec  1.27min
Code 
## First, note that we have the structured id = 'mrf1' and unstructured
## spatial effect id = 're2', also indicated in the model summary
summary(b_wheat_yield_MRF)

Call:
bamlss(formula = f_wheat_yield_MRF, family = "gaussian", data = Irrig_Rev_rice_wheat)
---
Family: gaussian 
Link function: mu = identity, sigma = log
*---
Formula mu:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + 
    g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", 
    xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.5740916 -1.1340813 -0.6630343  0.1609182     -1.359
variety_type_NMWV    0.3056829  0.2611116  0.3052558  0.3529448      0.303
g_q5305_irrig_times  0.4153445  0.3895350  0.4153230  0.4439696      0.411
nperha               0.0012336  0.0006668  0.0012288  0.0017957      0.001
p2o5perha            0.0031334  0.0020527  0.0031371  0.0042060      0.003
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                    Mean      2.5%       50%     97.5%
s(sowdate_fmt_wheat_day).tau21 4.037e-02 7.694e-05 3.583e-03 3.627e-01
s(sowdate_fmt_wheat_day).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(sowdate_fmt_wheat_day).edf   1.940e+00 1.015e+00 1.491e+00 4.943e+00
s(District,id='mrf1').tau21    3.294e-03 6.895e-05 9.061e-04 1.507e-02
s(District,id='mrf1').alpha    1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf      1.865e+01 3.933e+00 1.853e+01 3.206e+01
s(District,id='re2').tau21     4.793e+00 1.759e+00 4.533e+00 9.645e+00
s(District,id='re2').alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf       3.589e+01 3.577e+01 3.590e+01 3.596e+01
                               parameters
s(sowdate_fmt_wheat_day).tau21     15.401
s(sowdate_fmt_wheat_day).alpha         NA
s(sowdate_fmt_wheat_day).edf        8.309
s(District,id='mrf1').tau21         0.003
s(District,id='mrf1').alpha            NA
s(District,id='mrf1').edf          26.415
s(District,id='re2').tau21          6.085
s(District,id='re2').alpha             NA
s(District,id='re2').edf           35.897
---
Formula sigma:
---
sigma ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + 
    nperha + p2o5perha + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.7819694 -0.8953612 -0.7812289 -0.6696800     -0.316
variety_type_NMWV    0.2056530  0.1502403  0.2059671  0.2547450      0.198
g_q5305_irrig_times  0.0927109  0.0623800  0.0929975  0.1245421      0.098
nperha              -0.0011626 -0.0018926 -0.0011700 -0.0004521     -0.001
p2o5perha            0.0017864  0.0003726  0.0017712  0.0031635      0.002
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.8957 0.4855 0.9745     1
-
Smooth terms:
                                    Mean      2.5%       50%     97.5%
s(sowdate_fmt_wheat_day).tau21 8.472e-02 5.798e-04 4.162e-02 4.208e-01
s(sowdate_fmt_wheat_day).alpha 9.647e-01 7.699e-01 9.984e-01 1.000e+00
s(sowdate_fmt_wheat_day).edf   2.574e+00 1.067e+00 2.524e+00 4.535e+00
s(District,id='mrf1').tau21    6.406e-04 5.097e-05 3.717e-04 2.418e-03
s(District,id='mrf1').alpha    8.975e-01 4.521e-01 9.764e-01 1.000e+00
s(District,id='mrf1').edf      1.089e+01 2.084e+00 9.704e+00 2.304e+01
s(District,id='re2').tau21     1.986e-02 9.452e-03 1.872e-02 3.562e-02
s(District,id='re2').alpha     7.360e-01 1.839e-01 8.280e-01 1.000e+00
s(District,id='re2').edf       2.834e+01 2.429e+01 2.855e+01 3.129e+01
                               parameters
s(sowdate_fmt_wheat_day).tau21      0.035
s(sowdate_fmt_wheat_day).alpha         NA
s(sowdate_fmt_wheat_day).edf        2.415
s(District,id='mrf1').tau21         0.001
s(District,id='mrf1').alpha            NA
s(District,id='mrf1').edf          13.967
s(District,id='re2').tau21          0.203
s(District,id='re2').alpha             NA
s(District,id='re2').edf           34.551
---
Sampler summary:
-
DIC = 8219.792 logLik = -4070.916 pd = 77.9595
runtime = 76.61
---
Optimizer summary:
-
AICc = 8311.175 edf = 131.5535 logLik = -4020.075
logPost = -4029.064 nobs = 4537 runtime = 5.11
Code 
# Plot the nonlinear effect
plot(b_wheat_yield_MRF, model = "mu", term = "s(sowdate_fmt_wheat_day)")
Code 
## Now, to predict the spatial effects we set up new data.
nd <- data.frame("District" = unique(Irrig_Rev_rice_wheat$District))

## Predict for the structured spatial effects.
p_str_wheat_yield_MRF <- predict(b_wheat_yield_MRF, newdata = nd, term = "s(District,id='mrf1')", intercept = FALSE)

## And the unstructured spatial effect.
p_unstr_wheat_yield_MRF <- predict(b_wheat_yield_MRF, newdata = nd, term = "s(District,id='re2')", intercept = FALSE)

## MRF smooth effect.
plotmap(India_aoi_sp_bnd,
    x = p_str_wheat_yield_MRF$mu, id = nd$District,
    main = expression(mu), title = "Structured spatial effect"
)
Code 
plotmap(India_aoi_sp_bnd,
    x = p_str_wheat_yield_MRF$sigma, id = nd$District,
    main = expression(sigma), title = "Structured spatial effect"
)
Code 
## Random effects.
plotmap(India_aoi_sp_bnd,
    x = p_unstr_wheat_yield_MRF$mu, id = nd$District, title = "Unstructured spatial effect"
)
Code 
plotmap(India_aoi_sp_bnd,
    x = p_unstr_wheat_yield_MRF$sigma, id = nd$District, title = "Unstructured spatial effect"
)
Code 
library(sp)
library(mgcv)
library(bamlss)

# Multivariate geoadditive model
# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")
# install_github("https://github.com/meteosimon/mvnchol")
library(mvnchol)

# Remove NAs in the monsoon variables
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$onset_2017)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$monsoon_onset_dev)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$median_onset_82_15)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$sd_onset_82_15)))

# Rice first stage
f_sow_rice_1st_stage <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)

f_sow_rice_1st_stage_MRF <- bamlss(f_sow_rice_1st_stage, data = Irrig_Rev_rice_wheat, family = "gaussian")
AICc 38862.03 logPost -165306. logLik -19371.2 edf 59.010 eps 0.1226 iteration   1
AICc 34841.08 logPost -31072.7 logLik -17339.3 edf 79.788 eps 0.0980 iteration   2
AICc 33498.94 logPost -18399.0 logLik -16648.2 edf 98.970 eps 0.0632 iteration   3
AICc 33180.98 logPost -17093.4 logLik -16491.0 edf 97.287 eps 0.0297 iteration   4
AICc 33096.83 logPost -16986.5 logLik -16452.6 edf 93.751 eps 0.0069 iteration   5
AICc 33056.53 logPost -16935.2 logLik -16435.7 edf 90.674 eps 0.0017 iteration   6
AICc 33040.37 logPost -16899.0 logLik -16429.8 edf 88.587 eps 0.0009 iteration   7
AICc 33035.63 logPost -16869.0 logLik -16428.1 edf 87.868 eps 0.0003 iteration   8
AICc 33028.91 logPost -16845.9 logLik -16427.6 edf 85.110 eps 0.0001 iteration   9
AICc 33028.60 logPost -16840.2 logLik -16427.5 edf 85.077 eps 0.0000 iteration  10
AICc 33028.60 logPost -16840.2 logLik -16427.5 edf 85.077 eps 0.0000 iteration  10
elapsed time:  4.89sec
Starting the sampler...

|                    |   0%  1.39min
|*                   |   5%  1.36min  4.31sec
|**                  |  10%  1.18min  7.89sec
|***                 |  15%  1.15min 12.15sec
|****                |  20%  1.10min 16.47sec
|*****               |  25%  1.03min 20.58sec
|******              |  30% 57.38sec 24.59sec
|*******             |  35% 53.28sec 28.69sec
|********            |  40% 48.93sec 32.62sec
|*********           |  45% 43.19sec 35.34sec
|**********          |  50% 38.26sec 38.26sec
|***********         |  55% 33.59sec 41.06sec
|************        |  60% 30.39sec 45.58sec
|*************       |  65% 26.62sec 49.44sec
|**************      |  70% 22.64sec 52.83sec
|***************     |  75% 19.05sec 57.14sec
|****************    |  80% 15.28sec  1.02min
|*****************   |  85% 11.52sec  1.09min
|******************  |  90%  7.65sec  1.15min
|******************* |  95%  3.81sec  1.21min
|********************| 100%  0.00sec  1.26min
Code 
fitted_f_sow_rice_1st_stage_MRF <- f_sow_rice_1st_stage_MRF$fitted

Irrig_Rev_rice_wheat$Res_rice_sow <- Irrig_Rev_rice_wheat$sowdate_fmt_rice_day - fitted_f_sow_rice_1st_stage_MRF$mu

# Wheat first stage
f_sow_wheat_1st_stage <- list(
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(gw_2018) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)

f_sow_wheat_1st_stage_MRF <- bamlss(f_sow_wheat_1st_stage, data = Irrig_Rev_rice_wheat, family = "gaussian")
AICc 49188.65 logPost -488850. logLik -24543.4 edf 50.285 eps 0.2566 iteration   1
AICc 44824.42 logPost -64385.9 logLik -22386.0 edf 26.046 eps 0.0753 iteration   2
AICc 40409.05 logPost -24539.6 logLik -20136.0 edf 67.444 eps 0.0824 iteration   3
AICc 36648.46 logPost -18973.1 logLik -18246.6 edf 76.219 eps 0.0789 iteration   4
AICc 33900.77 logPost -17313.2 logLik -16871.9 edf 77.104 eps 0.0796 iteration   5
AICc 32636.01 logPost -16736.4 logLik -16238.1 edf 78.490 eps 0.0660 iteration   6
AICc 32419.72 logPost -16629.9 logLik -16129.2 edf 79.163 eps 0.0324 iteration   7
AICc 32413.31 logPost -16627.4 logLik -16125.7 edf 79.502 eps 0.0056 iteration   8
AICc 32412.92 logPost -16627.3 logLik -16125.4 edf 79.574 eps 0.0001 iteration   9
AICc 32412.79 logPost -16627.2 logLik -16125.3 edf 79.581 eps 0.0000 iteration  10
AICc 32412.79 logPost -16627.2 logLik -16125.3 edf 79.581 eps 0.0000 iteration  10
elapsed time:  2.08sec
Starting the sampler...

|                    |   0% 38.08sec
|*                   |   5% 54.72sec  2.88sec
|**                  |  10% 59.85sec  6.65sec
|***                 |  15% 53.95sec  9.52sec
|****                |  20% 49.20sec 12.30sec
|*****               |  25% 44.79sec 14.93sec
|******              |  30% 41.42sec 17.75sec
|*******             |  35% 39.24sec 21.13sec
|********            |  40% 36.22sec 24.15sec
|*********           |  45% 33.39sec 27.32sec
|**********          |  50% 31.19sec 31.19sec
|***********         |  55% 28.46sec 34.79sec
|************        |  60% 25.27sec 37.91sec
|*************       |  65% 21.87sec 40.61sec
|**************      |  70% 18.45sec 43.05sec
|***************     |  75% 15.01sec 45.04sec
|****************    |  80% 11.77sec 47.08sec
|*****************   |  85%  8.68sec 49.16sec
|******************  |  90%  5.70sec 51.27sec
|******************* |  95%  2.82sec 53.58sec
|********************| 100%  0.00sec 55.61sec
Code 
fitted_f_sow_wheat_1st_stage_MRF <- f_sow_wheat_1st_stage_MRF$fitted

Irrig_Rev_rice_wheat$Res_wheat_sow <- Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day - fitted_f_sow_wheat_1st_stage_MRF$mu


# Multivariate with sowing dates as endogenous
f_rice_wheat_yield_MRF_corr_endo <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lamdiag1 ~ 1,
    lamdiag2 ~ 1,
    lamdiag3 ~ 1,
    lamdiag4 ~ 1,
    lambda12 ~ 1,
    lambda13 ~ 1,
    lambda14 ~ 1,
    lambda23 ~ 1,
    lambda24 ~ s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lambda34 ~ 1
)

multivariate_geo_sow_MRF_corr_endo <- bamlss(f_rice_wheat_yield_MRF_corr_endo, type = "modified", family = mvnchol_bamlss(k = 4), data = Irrig_Rev_rice_wheat)
AICc 349476.1 logPost -2543606 logLik -174262. edf 430.14 eps 1.0000 iteration   1
AICc 157743.5 logPost -2447306 logLik -78382.3 edf 441.55 eps 0.1825 iteration   2
AICc 104307.5 logPost -2420173 logLik -51667.3 edf 439.10 eps 0.2070 iteration   3
AICc 89316.74 logPost -2412636 logLik -44189.9 edf 424.43 eps 0.1254 iteration   4
AICc 86452.65 logPost -2411216 logLik -42768.9 edf 415.34 eps 0.0628 iteration   5
AICc 86231.98 logPost -2411101 logLik -42661.9 edf 412.53 eps 0.0202 iteration   6
AICc 86215.73 logPost -2411097 logLik -42659.1 edf 408.13 eps 0.0027 iteration   7
AICc 86213.14 logPost -2411094 logLik -42658.2 edf 407.83 eps 0.0007 iteration   8
AICc 86211.64 logPost -2411092 logLik -42657.7 edf 407.62 eps 0.0003 iteration   9
AICc 86210.65 logPost -2411090 logLik -42657.4 edf 407.43 eps 0.0002 iteration  10
AICc 86205.44 logPost -2411208 logLik -42657.3 edf 405.37 eps 0.0001 iteration  11
AICc 86205.08 logPost -2411207 logLik -42657.3 edf 405.26 eps 0.0001 iteration  12
AICc 86202.45 logPost -2412410 logLik -42657.2 edf 404.20 eps 0.0001 iteration  13
AICc 86202.37 logPost -2412333 logLik -42657.2 edf 404.18 eps 0.0001 iteration  14
AICc 86202.37 logPost -2412333 logLik -42657.2 edf 404.18 eps 0.0001 iteration  14
elapsed time: 55.35sec
Starting the sampler...

|                    |   0%  8.03min
|*                   |   5%  7.75min 24.47sec
|**                  |  10%  7.15min 47.69sec
|***                 |  15%  6.83min  1.21min
|****                |  20%  6.57min  1.64min
|*****               |  25%  6.39min  2.13min
|******              |  30%  6.05min  2.59min
|*******             |  35%  5.64min  3.04min
|********            |  40%  5.27min  3.51min
|*********           |  45%  4.84min  3.96min
|**********          |  50%  4.42min  4.42min
|***********         |  55%  4.01min  4.90min
|************        |  60%  3.66min  5.49min
|*************       |  65%  3.21min  5.96min
|**************      |  70%  2.79min  6.51min
|***************     |  75%  2.35min  7.04min
|****************    |  80%  1.87min  7.50min
|*****************   |  85%  1.40min  7.95min
|******************  |  90% 56.05sec  8.41min
|******************* |  95% 28.00sec  8.87min
|********************| 100%  0.00sec  9.32min
Code 
summary(multivariate_geo_sow_MRF_corr_endo)

Call:
bamlss(formula = f_rice_wheat_yield_MRF_corr_endo, family = mvnchol_bamlss(k = 4), 
    data = Irrig_Rev_rice_wheat, type = "modified")
---
Family: mvnchol 
Link function: mu1 = identity, mu2 = identity, mu3 = identity, mu4 = identity, lamdiag1 = log, lamdiag2 = log, lamdiag3 = log, lamdiag4 = log, lambda12 = identity, lambda13 = identity, lambda14 = identity, lambda23 = identity, lambda24 = identity, lambda34 = identity
*---
Formula mu1:
---
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + 
    s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + 
    s(sd_onset_82_15) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              68.6313 65.0608 69.0601 70.4619      2.030
rice_duration_class_long -1.3633 -2.0104 -1.3634 -0.7134     -1.341
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(gw_2017).tau21            2.095e+00 1.260e-04 1.210e-02 2.125e+01    428.875
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              1.326e+00 9.987e-01 1.008e+00 3.388e+00      7.459
s(onset_2017).tau21         7.779e+00 1.169e-04 2.620e-01 7.322e+01   3370.678
s(onset_2017).alpha         1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(onset_2017).edf           1.790e+00 9.986e-01 1.275e+00 4.686e+00      8.245
s(monsoon_onset_dev).tau21  7.665e+00 2.474e-04 1.551e-01 7.089e+01     16.446
s(monsoon_onset_dev).alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(monsoon_onset_dev).edf    1.930e+00 9.996e-01 1.214e+00 5.147e+00      3.687
s(median_onset_82_15).tau21 2.959e+01 8.439e-04 1.024e+01 1.830e+02     18.779
s(median_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(median_onset_82_15).edf   3.182e+00 1.001e+00 3.276e+00 6.208e+00      3.766
s(sd_onset_82_15).tau21     1.316e+01 1.644e-03 1.563e+00 8.969e+01      0.537
s(sd_onset_82_15).alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(sd_onset_82_15).edf       2.170e+00 1.002e+00 1.801e+00 5.006e+00      1.379
s(District,id='mrf1').tau21 3.089e+01 1.575e+01 2.914e+01 5.586e+01   1484.368
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   3.430e+01 3.394e+01 3.431e+01 3.457e+01     34.979
s(District,id='re2').tau21  1.606e+04 1.019e+04 1.539e+04 2.517e+04      1.087
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.599e+01 3.599e+01 3.599e+01 3.600e+01     35.240
---
Formula mu2:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + 
    s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + 
    s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              2.16281 0.43171 1.97757 3.90909      0.000
rice_duration_class_long 0.10345 0.03601 0.10329 0.17645      0.096
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(Res_rice_sow).tau21             2.939e+00 7.670e-02 1.377e+00 1.492e+01
s(Res_rice_sow).alpha             1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(Res_rice_sow).edf               5.543e+00 2.630e+00 5.592e+00 8.087e+00
s(g_q5305_irrig_times_rice).tau21 1.715e+00 2.310e-01 1.227e+00 5.434e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf   6.171e+00 4.441e+00 6.195e+00 7.609e+00
s(nperha_rice).tau21              7.856e-01 5.357e-03 3.042e-01 3.674e+00
s(nperha_rice).alpha              1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf                3.921e+00 1.288e+00 3.786e+00 6.668e+00
s(p2o5perha_rice).tau21           2.061e-01 1.334e-04 3.188e-02 1.348e+00
s(p2o5perha_rice).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf             2.348e+00 1.010e+00 2.059e+00 5.205e+00
s(District,id='mrf1').tau21       3.043e-03 8.435e-05 1.067e-03 1.628e-02
s(District,id='mrf1').alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf         1.462e+01 2.134e+00 1.342e+01 3.017e+01
s(District,id='re2').tau21        4.974e+00 1.229e-01 4.138e+00 1.460e+01
s(District,id='re2').alpha        1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf          3.529e+01 3.323e+01 3.577e+01 3.593e+01
                                  parameters
s(Res_rice_sow).tau21                 16.244
s(Res_rice_sow).alpha                     NA
s(Res_rice_sow).edf                    8.170
s(g_q5305_irrig_times_rice).tau21      1.698
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        6.538
s(nperha_rice).tau21                   0.991
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     5.103
s(p2o5perha_rice).tau21                0.172
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  3.239
s(District,id='mrf1').tau21            0.035
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             32.259
s(District,id='re2').tau21             8.948
s(District,id='re2').alpha                NA
s(District,id='re2').edf              35.894
---
Formula mu3:
---
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + 
    s(District, bs = "mrf", xt = list(penalty = K)) + s(District, 
    bs = "re")
-
Parametric coefficients:
                     Mean    2.5%     50%   97.5% parameters
(Intercept)       114.655 109.483 114.481 120.587      3.830
variety_type_NMWV  -3.072  -3.742  -3.074  -2.460     -3.073
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(harvest_day_rice).tau21     228.435    44.388   160.237   704.149   1035.210
s(harvest_day_rice).alpha       1.000     1.000     1.000     1.000         NA
s(harvest_day_rice).edf         5.934     4.524     5.898     7.396      8.391
s(District,id='mrf1').tau21    54.003     9.448    45.568   154.095   4326.923
s(District,id='mrf1').alpha     1.000     1.000     1.000     1.000         NA
s(District,id='mrf1').edf      34.590    33.755    34.685    34.847     34.994
s(District,id='re2').tau21  49904.116 32148.970 48241.930 80332.511      1.087
s(District,id='re2').alpha      1.000     1.000     1.000     1.000         NA
s(District,id='re2').edf       35.999    35.998    35.999    35.999     35.240
---
Formula mu4:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + 
    g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", 
    xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.3451667 -0.9637935 -0.4133652  0.3762429     -1.419
variety_type_NMWV    0.3422512  0.2824899  0.3429785  0.3966410      0.337
g_q5305_irrig_times  0.4219196  0.3932712  0.4221273  0.4494278      0.424
nperha               0.0015496  0.0009547  0.0015550  0.0021134      0.002
p2o5perha            0.0025587  0.0014420  0.0025588  0.0036574      0.002
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(Res_wheat_sow).tau21      1.861e-02 7.837e-05 1.649e-03 1.452e-01      0.238
s(Res_wheat_sow).alpha      1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_wheat_sow).edf        1.607e+00 1.011e+00 1.202e+00 3.875e+00      4.795
s(District,id='mrf1').tau21 1.002e-02 5.295e-05 5.635e-03 4.258e-02      0.010
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   2.528e+01 2.963e+00 2.938e+01 3.356e+01     31.282
s(District,id='re2').tau21  3.750e+00 1.161e+00 3.640e+00 7.588e+00      5.082
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.586e+01 3.568e+01 3.588e+01 3.594e+01     35.915
---
Formula lamdiag1:
---
lamdiag1 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.216 -2.236 -2.216 -2.194     -2.207
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9896 0.9035 1.0000     1
---
Formula lamdiag2:
---
lamdiag2 ~ 1
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) 0.06651 0.04523 0.06601 0.08788      0.074
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9913 0.9196 1.0000     1
---
Formula lamdiag3:
---
lamdiag3 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.145 -2.167 -2.145 -2.126     -2.141
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9896 0.9163 1.0000     1
---
Formula lamdiag4:
---
lamdiag4 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) 0.5178 0.4983 0.5173 0.5394      0.527
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9910 0.9175 1.0000     1
---
Formula lambda12:
---
lambda12 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept) -0.011213 -0.045104 -0.007393  0.028365     -0.004
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda13:
---
lambda13 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.007309 0.003664 0.007292 0.010696      0.008
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda14:
---
lambda14 ~ 1
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) 0.01113 0.00202 0.01088 0.01917      0.011
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda23:
---
lambda23 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept)  0.023714 -0.007921  0.023630  0.055818      0.027
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda24:
---
lambda24 ~ s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) -0.2382 -0.2908 -0.2359 -0.1996     -0.036
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(District,id='mrf1').tau21 2.387e-03 1.004e-04 2.022e-03 6.647e-03      0.005
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   1.479e+01 1.863e+00 1.601e+01 2.430e+01     22.153
s(District,id='re2').tau21  9.553e-03 1.586e-04 6.116e-03 3.340e-02      0.024
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    1.374e+01 6.443e-01 1.415e+01 2.686e+01     24.454
---
Formula lambda34:
---
lambda34 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) -0.07590 -0.11822 -0.08111 -0.01156     -0.009
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Sampler summary:
-
DIC = 85815.46 logLik = -42797.02 pd = 221.4134
runtime = 561.25
---
Optimizer summary:
-
AICc = 86202.37 edf = 404.1805 logLik = -42657.27
logPost = -2412333 nobs = 4527 runtime = 55.35
Code 
nd <- data.frame("District" = unique(Irrig_Rev_rice_wheat$District))

# Focusing on the cross-equation correlations

## Predict for the structured spatial effects.
p_str_multivariate_geo_sow_MRF_corr_endo_rice_y <- predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term = "s(District,id='mrf1')", intercept = FALSE)

p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt <- as.data.frame(p_str_multivariate_geo_sow_MRF_corr_endo_rice_y)


## And the unstructured spatial effect.
p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y <- predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term = "s(District,id='re2')", intercept = FALSE)

# Rice sowing spatial equation

plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main = expression(mu(rice_sowing)), title = "Structured spatial effect")
Code 
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main = expression(mu(rice_sowing)), title = "Unstructured spatial effect")
Code 
# Rice yield spatial equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main = expression(mu(rice_yield)), title = "Structured spatial effect")
Code 
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main = expression(mu(rice_yield)), title = "Unstructured spatial effect")
Code 
# Wheat sowing equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main = expression(mu(wheat_sowing)), title = "Structured spatial effect")
Code 
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main = expression(mu(wheat_sowing)), title = "Unstructured spatial effect")
Code 
# Wheat yield spatial equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main = expression(mu(wheat_yield)), title = "Structured spatial effect")
Code 
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main = expression(mu(wheat_yield)), title = "Unstructured spatial effect")
Code 
# Focusing on the cross-equation correlations

## MRF smooth effect.
plotmap(India_aoi_sp_bnd,
    x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District,
    main = expression(lambda(rice, wheat)), title = "Structured spatial effect"
)
Code 
## Random effects.
plotmap(India_aoi_sp_bnd,
    x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District, main = expression(lambda(rice, wheat)), title = "Unstructured spatial effect"
)
Code 
# Rice sowing equation : Non-linear relationships
# s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15)
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(gw_2017)")
Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(onset_2017)")
Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(monsoon_onset_dev)")
Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(median_onset_82_15)")
Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(sd_onset_82_15)")
Code 
# Rice yield equation
# s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice)

plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(Res_rice_sow)")
Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(g_q5305_irrig_times_rice)")
Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(nperha_rice)")
Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(p2o5perha_rice)")
Code 
# Wheat sowing equation
# s(harvest_day_rice) + s(gw_2018)
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(harvest_day_rice)")
Code 
# (multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(gw_2018)")

# Wheat yield equation
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu4", term = "s(Res_wheat_sow)")
Code 
# Fitted values
multivariate_geo_sow_MRF_corr_endo_fitted_values <- multivariate_geo_sow_MRF_corr_endo$fitted

multivariate_geo_sow_MRF_corr_endo_fitted_values <- as.data.frame(multivariate_geo_sow_MRF_corr_endo_fitted_values)

# Merge the fitted results to the data and export



Irrig_Rev_rice_wheat_Mult_Res <- cbind(Irrig_Rev_rice_wheat, multivariate_geo_sow_MRF_corr_endo_fitted_values)


summary(lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res))

Call:
lm(formula = sowdate_fmt_rice_day ~ mu1, data = Irrig_Rev_rice_wheat_Mult_Res)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.574  -5.397  -0.269   4.986  44.985 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.03383    3.41385   -0.01    0.992    
mu1          1.00018    0.01771   56.48   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.093 on 4525 degrees of freedom
Multiple R-squared:  0.4135,    Adjusted R-squared:  0.4134 
F-statistic:  3190 on 1 and 4525 DF,  p-value: < 2.2e-16
Code 
summary(lm(b_grain_yield_ton_per_ha_rice ~ l_ton_per_hectare, Irrig_Rev_rice_wheat_Mult_Res))

Call:
lm(formula = b_grain_yield_ton_per_ha_rice ~ l_ton_per_hectare, 
    data = Irrig_Rev_rice_wheat_Mult_Res)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2891 -0.7976 -0.0052  0.7454  3.3894 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        3.34109    0.05974   55.93   <2e-16 ***
l_ton_per_hectare  0.22766    0.01958   11.63   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.08 on 4525 degrees of freedom
Multiple R-squared:  0.02901,   Adjusted R-squared:  0.0288 
F-statistic: 135.2 on 1 and 4525 DF,  p-value: < 2.2e-16
Code 
plot((lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res)))
Code 
Irrig_Rev_rice_wheat_Mult_Res_sp <- Irrig_Rev_rice_wheat_Mult_Res
coordinates(Irrig_Rev_rice_wheat_Mult_Res_sp) <- c("o_largest_plot_gps_longitude", "o_largest_plot_gps_latitude")

# Map the correlations
# library(modelsummary)
# Irrig_Rev_rice_wheat_Mult_Res_dist <- datasummary(Heading("District") * District ~ Heading("N obs") * N + Heading("%") * Percent() + lambda24 * (Mean + SD), data = Irrig_Rev_rice_wheat_Mult_Res, output = "data.frame")

# library(dplyr)
# Irrig_Rev_rice_wheat_Mult_Res_dist <- rename(Irrig_Rev_rice_wheat_Mult_Res_dist, "Mean_Rice_Wheat_Rho" = "Mean")
# Irrig_Rev_rice_wheat_Mult_Res_dist <- rename(Irrig_Rev_rice_wheat_Mult_Res_dist, "SD_Rice_Wheat_Rho" = "SD")

# Irrig_Rev_rice_wheat_Mult_Res_dist <- subset(Irrig_Rev_rice_wheat_Mult_Res_dist, Irrig_Rev_rice_wheat_Mult_Res_dist$District != "Purnia")

# rice_wheat_yield_rho_dist_sf <- merge(India_aoi_sf, Irrig_Rev_rice_wheat_Mult_Res_dist, by = "District")

# rice_wheat_yield_rho_dist_sf$Mean_Rice_Wheat_Rho <- as.numeric(rice_wheat_yield_rho_dist_sf$Mean_Rice_Wheat_Rho)
# library(mapview)
# mapview(rice_wheat_yield_rho_dist_sf, zcol = "Mean_Rice_Wheat_Rho", layer.name = "Rice wheat equation correlation")
# library(sf)
# rice_wheat_yield_rho_dist_sf_sp <- as_Spatial(rice_wheat_yield_rho_dist_sf)
# library(tmap)
# tmap_mode("view")
# rice_wheat_yield_rho_dist_sf_sp_map <- tm_shape(rice_wheat_yield_rho_dist_sf_sp) +
#     tm_polygons(col = "Mean_Rice_Wheat_Rho", title = "Rice wheat equation correlation", style = "quantile") +
#     tm_layout(legend.outside = TRUE)

# tmap_save(rice_wheat_yield_rho_dist_sf_sp_map, "figures/rice_wheat_yield_rho_dist_sf_sp_map .png")
Code 
f_rice_wheat_yield_MRF_corr_endo_General <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lamdiag1 ~ 1,
    lamdiag2 ~ 1,
    lamdiag3 ~ 1,
    lamdiag4 ~ 1,
    lambda12 ~ 1,
    lambda13 ~ 1,
    lambda14 ~ 1,
    lambda23 ~ 1,
    lambda24 ~ s(gw_2017)+s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lambda34 ~ 1
)

multivariate_geo_sow_MRF_corr_endo_General <- bamlss(f_rice_wheat_yield_MRF_corr_endo_General, type = "modified", family = mvnchol_bamlss(k = 4), data = Irrig_Rev_rice_wheat)
AICc 349479.0 logPost -2543604 logLik -174260. edf 432.85 eps 1.0000 iteration   1
AICc 157745.3 logPost -2447296 logLik -78381.2 edf 443.21 eps 0.1854 iteration   2
AICc 104308.9 logPost -2420161 logLik -51666.0 edf 440.78 eps 0.2158 iteration   3
AICc 89318.46 logPost -2412624 logLik -44188.5 edf 426.28 eps 0.1451 iteration   4
AICc 86454.46 logPost -2411205 logLik -42767.5 edf 417.25 eps 0.0645 iteration   5
AICc 86233.78 logPost -2411090 logLik -42660.5 edf 414.43 eps 0.0206 iteration   6
AICc 86217.50 logPost -2411085 logLik -42657.7 edf 410.02 eps 0.0028 iteration   7
AICc 86214.90 logPost -2411082 logLik -42656.8 edf 409.71 eps 0.0010 iteration   8
AICc 86213.39 logPost -2411080 logLik -42656.3 edf 409.50 eps 0.0007 iteration   9
AICc 86212.40 logPost -2411078 logLik -42656.0 edf 409.31 eps 0.0045 iteration  10
AICc 86207.19 logPost -2411196 logLik -42655.9 edf 407.25 eps 0.0012 iteration  11
AICc 86206.82 logPost -2411195 logLik -42655.9 edf 407.14 eps 0.0003 iteration  12
AICc 86204.18 logPost -2412413 logLik -42655.8 edf 406.07 eps 0.0003 iteration  13
AICc 86204.09 logPost -2412336 logLik -42655.8 edf 406.05 eps 0.0007 iteration  14
AICc 86204.03 logPost -2412276 logLik -42655.8 edf 406.03 eps 0.0002 iteration  15
AICc 86203.96 logPost -2412276 logLik -42655.8 edf 406.01 eps 0.0001 iteration  16
AICc 86203.91 logPost -2412276 logLik -42655.8 edf 405.99 eps 0.0001 iteration  17
AICc 86203.91 logPost -2412276 logLik -42655.8 edf 405.99 eps 0.0001 iteration  17
elapsed time:  1.07min
Starting the sampler...

|                    |   0%  6.31min
|*                   |   5%  5.89min 18.61sec
|**                  |  10%  5.66min 37.71sec
|***                 |  15%  5.40min 57.21sec
|****                |  20%  5.20min  1.30min
|*****               |  25%  5.01min  1.67min
|******              |  30%  4.78min  2.05min
|*******             |  35%  4.49min  2.42min
|********            |  40%  4.18min  2.79min
|*********           |  45%  3.86min  3.16min
|**********          |  50%  3.53min  3.53min
|***********         |  55%  3.19min  3.90min
|************        |  60%  2.85min  4.27min
|*************       |  65%  2.50min  4.64min
|**************      |  70%  2.14min  5.00min
|***************     |  75%  1.79min  5.38min
|****************    |  80%  1.44min  5.74min
|*****************   |  85%  1.08min  6.13min
|******************  |  90% 43.32sec  6.50min
|******************* |  95% 21.68sec  6.87min
|********************| 100%  0.00sec  7.24min
Code 
summary(multivariate_geo_sow_MRF_corr_endo_General)

Call:
bamlss(formula = f_rice_wheat_yield_MRF_corr_endo_General, family = mvnchol_bamlss(k = 4), 
    data = Irrig_Rev_rice_wheat, type = "modified")
---
Family: mvnchol 
Link function: mu1 = identity, mu2 = identity, mu3 = identity, mu4 = identity, lamdiag1 = log, lamdiag2 = log, lamdiag3 = log, lamdiag4 = log, lambda12 = identity, lambda13 = identity, lambda14 = identity, lambda23 = identity, lambda24 = identity, lambda34 = identity
*---
Formula mu1:
---
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + 
    s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + 
    s(sd_onset_82_15) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              71.3421 66.6505 71.2306 76.3630      2.031
rice_duration_class_long -1.3725 -2.0659 -1.3734 -0.7002     -1.342
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(gw_2017).tau21            2.783e+01 3.293e-02 6.432e+00 1.946e+02    428.733
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              2.672e+00 1.023e+00 2.494e+00 5.730e+00      7.459
s(onset_2017).tau21         1.876e+00 6.997e-05 1.240e-02 2.299e+01   3363.637
s(onset_2017).alpha         1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(onset_2017).edf           1.333e+00 9.976e-01 1.016e+00 3.611e+00      8.244
s(monsoon_onset_dev).tau21  7.569e+00 7.212e-04 4.383e-01 5.875e+01     15.486
s(monsoon_onset_dev).alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(monsoon_onset_dev).edf    2.046e+00 1.001e+00 1.480e+00 5.042e+00      3.633
s(median_onset_82_15).tau21 1.339e+01 1.459e-04 4.668e-02 1.026e+02     18.784
s(median_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(median_onset_82_15).edf   2.045e+00 9.990e-01 1.063e+00 5.504e+00      3.766
s(sd_onset_82_15).tau21     3.826e+00 1.450e-04 2.333e-02 4.081e+01      0.554
s(sd_onset_82_15).alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(sd_onset_82_15).edf       1.480e+00 9.990e-01 1.023e+00 4.173e+00      1.388
s(District,id='mrf1').tau21 7.390e+01 3.864e+01 7.183e+01 1.273e+02   1481.527
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   3.464e+01 3.443e+01 3.465e+01 3.479e+01     34.979
s(District,id='re2').tau21  1.519e+04 9.677e+03 1.455e+04 2.354e+04      1.087
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.599e+01 3.599e+01 3.599e+01 3.600e+01     35.240
---
Formula mu2:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + 
    s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + 
    s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              1.33448 0.07265 1.08753 3.51466      0.001
rice_duration_class_long 0.09706 0.03197 0.09722 0.15912      0.095
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(Res_rice_sow).tau21             2.419e+00 9.364e-02 1.038e+00 1.288e+01
s(Res_rice_sow).alpha             1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(Res_rice_sow).edf               5.310e+00 2.798e+00 5.241e+00 7.992e+00
s(g_q5305_irrig_times_rice).tau21 1.857e+00 2.474e-01 1.345e+00 6.468e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf   6.226e+00 4.500e+00 6.274e+00 7.728e+00
s(nperha_rice).tau21              7.115e-01 1.034e-02 3.575e-01 3.765e+00
s(nperha_rice).alpha              1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf                3.949e+00 1.483e+00 3.963e+00 6.671e+00
s(p2o5perha_rice).tau21           3.420e-01 1.193e-04 3.836e-02 2.565e+00
s(p2o5perha_rice).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf             2.583e+00 1.008e+00 2.155e+00 5.870e+00
s(District,id='mrf1').tau21       7.992e-03 8.542e-05 2.262e-03 4.889e-02
s(District,id='mrf1').alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf         1.855e+01 2.203e+00 1.895e+01 3.276e+01
s(District,id='re2').tau21        8.530e+00 4.015e-01 8.099e+00 2.101e+01
s(District,id='re2').alpha        1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf          3.576e+01 3.472e+01 3.588e+01 3.595e+01
                                  parameters
s(Res_rice_sow).tau21                 15.978
s(Res_rice_sow).alpha                     NA
s(Res_rice_sow).edf                    8.159
s(g_q5305_irrig_times_rice).tau21      1.718
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        6.551
s(nperha_rice).tau21                   0.985
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     5.097
s(p2o5perha_rice).tau21                0.172
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  3.240
s(District,id='mrf1').tau21            0.035
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             32.263
s(District,id='re2').tau21             8.982
s(District,id='re2').alpha                NA
s(District,id='re2').edf              35.894
---
Formula mu3:
---
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + 
    s(District, bs = "mrf", xt = list(penalty = K)) + s(District, 
    bs = "re")
-
Parametric coefficients:
                     Mean    2.5%     50%   97.5% parameters
(Intercept)       107.729 105.228 107.925 109.056      3.830
variety_type_NMWV  -3.071  -3.671  -3.072  -2.456     -3.073
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(harvest_day_rice).tau21     249.145    45.941   169.557   943.694   1033.860
s(harvest_day_rice).alpha       1.000     1.000     1.000     1.000         NA
s(harvest_day_rice).edf         5.867     4.288     5.819     7.667      8.391
s(District,id='mrf1').tau21   273.044   152.159   257.988   484.946   4322.922
s(District,id='mrf1').alpha     1.000     1.000     1.000     1.000         NA
s(District,id='mrf1').edf      34.913    34.840    34.914    34.977     34.994
s(District,id='re2').tau21  52391.176 32990.950 50597.124 79970.994      1.087
s(District,id='re2').alpha      1.000     1.000     1.000     1.000         NA
s(District,id='re2').edf       35.998    35.997    35.998    35.999     35.240
---
Formula mu4:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + 
    g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", 
    xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.7126342 -1.3986803 -0.8445395  0.2712965     -1.419
variety_type_NMWV    0.3424046  0.2887866  0.3427579  0.3931747      0.337
g_q5305_irrig_times  0.4231647  0.3972220  0.4229748  0.4504678      0.424
nperha               0.0015341  0.0009313  0.0015484  0.0021222      0.002
p2o5perha            0.0025301  0.0013402  0.0025688  0.0036608      0.003
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(Res_wheat_sow).tau21      2.950e-02 6.382e-05 2.565e-03 2.716e-01      0.233
s(Res_wheat_sow).alpha      1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_wheat_sow).edf        1.740e+00 1.009e+00 1.296e+00 4.531e+00      4.773
s(District,id='mrf1').tau21 2.583e-03 6.088e-05 5.571e-04 2.028e-02      0.010
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   1.662e+01 3.360e+00 1.473e+01 3.267e+01     31.294
s(District,id='re2').tau21  5.518e+00 1.334e+00 5.472e+00 1.160e+01      5.107
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.589e+01 3.571e+01 3.592e+01 3.596e+01     35.915
---
Formula lamdiag1:
---
lamdiag1 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.216 -2.237 -2.216 -2.196     -2.207
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9909 0.9239 1.0000     1
---
Formula lamdiag2:
---
lamdiag2 ~ 1
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) 0.06601 0.04535 0.06603 0.08673      0.074
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9907 0.9241 1.0000     1
---
Formula lamdiag3:
---
lamdiag3 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.145 -2.166 -2.145 -2.125     -2.141
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9905 0.9188 1.0000     1
---
Formula lamdiag4:
---
lamdiag4 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) 0.5179 0.4987 0.5176 0.5392      0.527
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9884 0.9012 1.0000     1
---
Formula lambda12:
---
lambda12 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept) -0.019167 -0.039310 -0.018658  0.002473     -0.004
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda13:
---
lambda13 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.006978 0.003265 0.007021 0.010635      0.008
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda14:
---
lambda14 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.013333 0.008215 0.013348 0.019123      0.011
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda23:
---
lambda23 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept)  0.023332 -0.008003  0.023810  0.053160      0.027
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda24:
---
lambda24 ~ s(gw_2017) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) -0.2375 -0.2912 -0.2359 -0.2011     -0.037
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(gw_2017).tau21            3.959e-02 8.403e-05 6.952e-03 2.810e-01      0.019
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              1.600e+00 1.005e+00 1.309e+00 3.396e+00      1.644
s(District,id='mrf1').tau21 3.633e-03 2.214e-04 3.565e-03 8.218e-03      0.005
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   1.832e+01 3.675e+00 1.986e+01 2.537e+01     22.402
s(District,id='re2').tau21  4.618e-03 1.066e-04 1.174e-03 2.527e-02      0.024
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    7.958e+00 4.446e-01 4.264e+00 2.514e+01     24.433
---
Formula lambda34:
---
lambda34 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) -0.04277 -0.16725 -0.04010  0.03504     -0.009
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Sampler summary:
-
DIC = 85820.91 logLik = -42798.93 pd = 223.049
runtime = 435.95
---
Optimizer summary:
-
AICc = 86203.92 edf = 405.9967 logLik = -42655.86
logPost = -2412277 nobs = 4527 runtime = 64.48
Code 
library(distreg.vis)

if (interactive()) {
    vis()
}